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Abstract 

A survey of linearized cosmological fluid equations with a number of different matter components 
is made. To begin with, the one-component case is reconsidered to illustrate some important 
mathematical and physical points rarely discussed in the literature. The work of some previous 
studies of two-component systems are examined and re-analyzed to point out some deficiencies of 
■ solutions, and further solutions and physical interpretation are then presented. This leads into 

a general two-component model with variable velocity dispersion parameters and mass density 
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CS| ' fractions of each component. The equations, applicable to both hot dark matter (HDM) and 
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•/^ ' cold dark matter (CDM) universes are solved in the long wavelength limit. This region is of 

. interest, because some modes in this range of wavenumbers are Jeans unstable. The mixture Jeans 

' wavenumber of the two-component system is introduced and interpreted, and the solutions are 

Qh' discussed, particularly in comparison to analogous solutions previously derived for plasma modes. 
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', structure formation is thought to have occurred. 
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I. INTRODUCTION 



The theory of structure formation in the Universe has become one of the most popu- 
lar and intensely studied topics in modern cosmology. Throughout the past century there 
has been an accumulating volume of work on the analytic investigation of the cosmologi- 
cal structure formation equations. The various approaches include both fluid and kinetic 
theory formulations. They principally consider the gravitational interaction of components 
of the cosmological medium, though sometimes other forms of interaction such as magnetic 
fields are also included (for some standard examples see e.g. . The analysis of these 

equations has employed ever more diverse and complicated techniques and approximation 
schemes to model increasingly realistic physical situations. This has been comprehensively 
supported and now superseded by large N-body simulations. The algorithms which govern 
these large numerical studies have grown progressively more refined and subtle, and are 
now producing very accurate and realistic results, which can be directly compared with 
observations (e.g.jbl). 

Despite the current trends in modern cosmological structure formation theory, much can 
still be learned from relatively simple analytic models. We consider such models, in the face 
of modern computing power, to analyze at a fundamental level some of the basic physical 
processes which cause the clustering observed in the Universe. This helps to isolate physical 
mechanisms difficult to discern numerically. In this paper our interest will focus on the 
linearized cosmological fluid equations. These equations have been used to build up the 
components of the cosmological density perturbation power spectrum, and must be evolved 
through the various stages of cosmological evolution, and over a large range of physical 
scales. There are several reasons for taking such an approach. The equations may be solved 
numerically to give detailed power spectra for the various cosmological models currently 
viable. The power spectra may then be used as initial data to evolve the large N-body 
simulations, which are ultimately compared to observations. The equations may also be 
used to build up a semi- quantitative picture of the evolution of the power spectrum. This 
can show how the various sized perturbations scale with respect to the Friedmann expansion 
parameter a during different epochs of the Universe, and give a direct insight into some of 
the fundamental physical processes operating to produce structure in the Universe (see e.g. 
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The evolution of perturbation modes with wavelength greater than the Hubble radius 
may be studied through a relativistic formulation of the perturbation equations, whereas for 
modes with wavelengths smaller than this radius, a Newtonian formalism suffices. There 
is a range of physical parameters and variety of differential equations describing the evo- 
lution processes of density fluctuations in the early Universe. This involves such elements 
as equations of state for energy constituents, and specification of the expansion parameter 
a by the Friedmann cosmological equations. As a consequence, there is a large wealth of 
literature on this subject, and most of the currently important techniques and results have 
been collected in some well known textbooks [l|, |2, y]. Some of the relatively complicated 
systems of equations have also been studied in the literature, and it is our goal to both 
review many of these studies, and to extend them in new directions, achieving some new 
unique results. 

In this paper we will only study the Newtonian limit of the linearized cosmological per- 
turbation equations, valid for density fiuctuations on scales well within the Hubble radius. 
Our main concern is with some mathematically more complicated multi-component models, 
which although not usually considered in standard power spectrum analysis, have realistic 
and interesting physical meaning. The primary concept of the Jeans gravitational instability 
Isll has been investigated in a static universe for multi-component models, to reveal the more 

nn 

complicated structure of modes possible ja, l2| • This provides some interesting qualitative 
ideas about the possible mechanisms for structure formation, but the lack of an expanding 
background spacetime in the models leads to unrealistic solutions, exponential in form. The 
inclusion of cosmological expansion in the equations leads to the more realistic power law 
and logarithmic solutions, familiar from the standard power spectrum analysis. Previous 
work in this area has focused both on some particular models jsl, and on a more general 
classification of the equations and solutions for a ranee of parameter values (some of them 
o„„ of n.tHen,a«cal pK..« context. Analytic solution, to. .o.e 

of the most general cases of the equations considered above, which often have significant 
physical interest, have not been achieved. It is our aim here to rectify this situation and 
investigate a system of equations modeling a two-component fiuid in the matter dominated 
post-recombination era of an Einstein-deSitter universe. One of these components consists 
of baryons, and the other some form of nonrelativistic dark matter particles. Through this 
work we will amend what appear to be some errors in the previous general studies of Haubold 
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and Mathai 11 1. 

This paper makes a comprehensive study of the hnear perturbation equations for cos- 
mological fluids with gravitational instabilities with application to large scale structure for- 
mation. For a historical perspective we note that Lifshitz concluded that gravitational 
instability could not be responsible for the formation of structure in the Universe. The 
correct conclusion, that gravitational instability suffices, was pointed out by Novikov jisl ]. 
As our work details the mathematical structure of the appropriate equations describing cos- 
mological structure formation, we note the nice series of papers by Ratra and Peebles Q| 
directed to understanding the applications of special functions to the problem of gravita- 
iional instability in cosmological models. We further note the nice series by Buchert et al 
concerning analytic results and their relevance to observational cosmology. 

We will also make a comparison with work done in cosmological plasma physics in an 



Einstein-deSitter background 
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j . This is interesting due to the mathematically very 
similar form of fluid equations for both type of systems. This similarity is largely due to 
the similarity of the electromagnetic and gravitational forces. In this paper we will analyze 
the long wavelength region of the solutions. This corresponds to the Jeans unstable region 
of parameter space, and requires use of Frobenius methods of expansion of the differential 
equations. In a follow-up paper jl^ we will investigate the short wavelength region of the 
solutions, which will require a WKB approximation scheme to be developed. 

In these papers we take explictly the temperature relationship T ~ l/a(t)^, where a{t) is 
the radius of the Universe. We give here the explanation why this is so. Following standard 
textbook material in Padmanabhan [l| (as summarized in equation (3.118)) and Peebles 
(1993, page 179) we find that the baryons follow this relationship when 

l + z = 142(nb/iV0.024)2/^ 

Here has been scaled to the WMAP best fit value. Thus for redshifts below z ~ 140 the 
baryon temperature drops as 1/a? down to 2; ~ 20, where the Universe reionizes (probably 
in a patchy fashion) and the 1/a? scaling no longer holds. This 20 < 2; < 140 redshift region 
is important, as it is where early large scale structure formation is thought to have occurred. 

The paper is to be organized as follows. In Section 2 we introduce the most general 
cosmological density perturbation equations in the Newtonian approximation. We review 
and classify previous work on these equations to put our current work into context, showing 
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what has been achieved, what needs amendment, and where we will seek to expand current 
knowledge. In Section 3 we revert to the one-component equations, to illustrate some of 
the basic principles which will be important later in our analysis, and to reveal some ap- 
parently new results. This will enable us to begin to tackle the two-component problem in 
Section 4. In this section the CDM two-component model in an expanding universe will 
be investigated. This has ties with the previous work cited, and we will demonstrate the 
limitations of the existing formalism here. We present new results apparently overlooked in 
the work of Haubold and Mathai After this, we are ready to study the most general 
baryonic and dark matter equations in Section 5, where we will consider the long wavelength 
approximation, applicable to either HDM or CDM. This is followed by our conclusions in 
Section 6. 



II. A CLASSIFICATION OF COSMOLOGICAL DENSITY PERTURBATION 
EQUATIONS 

As discussed in Section 1, there are a vast spectrum of equations describing cosmological 
density perturbations in different physical regimes. We begin directly with the linearized 
Newtonian approach. The equations for an n-component system of nonrelativistic species is 
derived in all the standard texts. Given a density perturbation 5i in the z-th component of 
the mass density p,: 

5,(r,t) = ^, (2.1) 

Pi 

it may be decomposed into its Fourier plane wave modes with wave vector k 

^i(r, t) = j SUt) exp(-ik • r)d?r. (2.2) 

Here r is the physical spatial coordinate, and t is cosmic time. 

To be able to solve the equations, the implicit time dependence of the physical variables 
needs to be removed. We will adopt the convention that barred variables will denote co- 
moving quantities, independent of time. Thus we define the comoving wave number k = ok. 
Using the Eulerian equations of motion describing a perfect fluid, a set of coupled second or- 
der equations for the Fourier modes 5i{t) (where we now drop the subscript k) are achieved: 

i=l 



Overdots will denote derivatives with respect to t. The above equations contain the sound 
velocity 

..=-ocp7 . (2.4) 

The sound velocity depends on the equation of state of the medium, and is in general 
dependent on time. We have introduced the specific heat ratio 7^, and assumed an equation 
of state of the form pi oc p]' . 

The introduction of a sound velocity implicitly assumes that the fluids involved are coUi- 
sional. This means that there are considerable interactions between the particles comprising 
each matter component. It is in fact generally assumed that dark matter is coUisionless, in 
which case a fluid equation is not strictly correct. The dark matter component would better 
be modeled by a distribution function satisfying the Vlasov equation. "Fluid-like" equations 
can still be derived in this case by taking velocity moments of the Vlasov equation, and iden- 
tifying the velocity dispersion with the above parameter fj. Thus although in the present 
paper we will refer to "sound velocities", this should be taken as a generic expression for a 
velocity dispersion parameter. Such an approach should work fine for CDM, but may neglect 
an important damping term found in HDM models. The complete analogous equation to 
fl2.3|) for HDM gives a fluid equation perspective on free-streaming, the phenomenon found 
in HDM models of neutrino-like matter. An approximate equation has been derived for 
a hot neutrino-like component by Setayeshgar [20[, which tends to wipe out perturbations 
below a certain scale (see also the lecture notes by Bertschinger [2^). In this work the exact 
Vlasov equation kinetic treatment was considered, and the usual Fermi-Dirac distribution 
function was replaced by a carefully chosen approximate form, allowing the conversion of 
the integro-differential equation into the following pure differential equation: 

— + —if 5, + -^5, = A-nG 2^ pA. (2.5) 

\ (X (XI (X 

^ ' 1=1 

Here the damping term 2kvi, gives rise to non-oscillating solutions heavily damped at short 
wavelengths. Thus an equation of the form ()2.3|) is not correct for HDM of a neutrino-like 
nature. We examine the general equations for both CDM and HDM, without specifying 
too carefully the exact nature of the dark component involved. This allows comparison of 
the results in this paper with previous work in the literature, which has also neglected this 
point. If some aspects of HDM models are poorly described by ()2.3j) . the equations are still 
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applicable to other two- component cosmological systems such as a hydrogen-helium gas not 
in equilibrium, where the lighter hydrogen component has a greater sound speed. 

At present (j2.3|) has been displayed in a quite general form, with an unspecified scale 
factor a, given by the Friedmann cosmological equation 
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General parameters describing the nature of the Universe in this equation are = 0, ±1, 
the spatial curvature, and A, the cosmological constant. It is difficult to make much progress 
without first becoming more specific about the energy content of the Universe. All the studies 



cited previously 



111 ] have only examined the Einstein-deSitter, matter dominated 
case, with various physical components and equations of state possible. The studies jlO. 



all of which are equivalent, make the pretense to study the radiation dominated era as well, 
but this is incorrect for the equations presented. It was explicitly assumed that a oc t^/^ 
in the scaling of the energy densities pi = Vti/ {Q-KGt^) and the velocity parameters, yet an 
allowance was made for a general Hubble expansion H = rjt^^. The general parameter t] can 
only be equal to | for the equations presented to be physically correct. The fiuid equations 
were formulated to allow for general equations of state, by writing the sound velocities such 
that both their magnitudes and time dependences were freely parameterized. A range of 
solutions were obtained for different cases of the parameters, and were generally classified 



by Meijer G-functions |22L |23| . We will show that the solutions found for a CDM and 
baryon model have been evaluated incorrectly, and will proceed to find their general exact 
representation. We will also proceed to investigate a more general dark matter and baryon 
problem than considered in any of the above. In j^, several subcases of the above mentioned 
studies were considered in some detail, and given a range of physical interpretations. The 
solutions were of a mathematically simpler nature, involving either Bessel functions or simple 
power law behavior. The investigations in Q concentrated on a three-component medium, 
involving baryons, CDM and photons. They incorrectly used the nonrelativistic Newtonian 
cosmological equations to model the photon component, so that the solutions, expressed in 
terms of Meijer G-functions, cannot be considered as physically relevant. 

Let us now make the choice of the matter dominated era of cosmological evolution in which 
to set ()2.3|1 . and in particular the post-recombination era, where baryons had decoupled from 
photons. This allows us to determine how the energy density and sound velocity scale with 
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respect to a, and consequently exhibit all explicit time dependences in the equations. We 
therefore introduce the comoving total background density po = a^po, and the constant 
ej = Pil Po, the fraction of the total mass density contributed by species i. This is distinct 
from the in general time dependent quantity VLi{t) = Pi/ Pc-, where pc is the critical density 
of the Universe 

We will consider a two-component fluid comprised of baryons (subscripted by B) and dark 
matter (subscripted by D). In the post-recombination era, the adiabatic speed of sound of 
species i assumes the following behavior: 

(xTi(x a-^ (2.8) 

where Tj is the temperature of the component. This prompts us to define the time indepen- 
dent quantity vf = a^vf. With these definitions, the linearized cosmological perturbation 
equations may be written as 

^ + + = ^{es5s + enSn). (2.10) 

The equations currently still represent a fairly general cosmological setting. The curvature 
parameter kc and cosmological constant A have not been specified, and control the behavior 
of a through the Friedmann equation ()2.6p . To see how these influence the evolution of the 
density perturbations, we transform the dependent variable from t to a. We also use (j2.6|) 
and another cosmological dynamics equation for the acceleration of a: 

a=lvrG^ + -a. (2.11) 
3 a2 3 ^ ^ 

This equation is derived in conjunction with the Friedmann equation by taking the spatial 
components of the Einstein equation. The cosmological perturbation equations are now 
able to be written in a form purely dependent on a, and parameterized explicitly by the 
cosmological dynamical constants: 

^gTrGpo + f «^ — kca)5"B + a^^^AixGpQ + ka? — 2kca)S'B 



'-^Sb - ^{eB5B + eoSo) = 0, (2.12) 
\7rGp0 + ^a^ — kca)6'^ + a~^{4TTGpo + Aa^ — 2kca)6'^ 

+ '-^Sd - ^{eB6B + eoSn) = 0. (2.13) 



In the above, a prime denotes differentiation with respect to a. 

A general analysis of ()2.12|) and ()2.13|) has not been attempted previously. To begin with, 
we can test for exactness of the equations (see j^J], pp. 92, 93), to determine whether a first 
integral exists. In general, a second order ordinary differential equation of the form 

AQ{x)y" + A^{x)y + A2{x)y = Q (2.14) 

[arbitrary functions Aj(x)] is exact if 

Ao-A; + A2 = 0. (2.15) 

Let us first consider the one-component example. This is the uncoupled case of ()2.12|) 
(e/) = 0, eB = 1), which gives 

K-X, + A, = ^^. (2.16) 

Thus the one-component equation is only exact for a pressureless gas vb = 0. This means 
that no closed form solution is possible, and approximations need to be made. We note 
that the pressureless one-component case has been studied extensively (e.g. Qj) for various 
values of the parameters. 

To make progress with the perturbation equations, and also to make contact with previous 
work in the literature, we need to make some assumptions about kc and A. The fee 7^ 
cases tend to be more complicated mathematically, as generally only parametric solutions 
can be found, where a is represented by hyperbolic functions {k^ = —1 open universe) or 
trigonometric functions {kc = 1 closed universe). Current observations, and the weight 
of theoretical tendencies in cosmology (e.g. = 1 as demanded by inflation) make the 
choice of flat universe kc = seem the most favorable. Q contains a contribution from A 
as well as matter components. The large amount of observational data now being analyzed, 
increasingly points to the existence of a cosmological constant comprising a major fraction of 
the energy density (see e.g. 25], [s^]), with a value of Q/<^ ^ 0.7. The Einstein-deSitter 



(A = 0, i7 = 1) model is generally not the model of choice anymore for detailed numerical 
studies in cosmology, however we do not make a claim that the solutions presented here 
are of an exact quantitative nature. Many other factors must also be taken into account 
when attempting to build up an exact, numerical model of structure formation. We wish 
to correct and extend some previous results, as well as perform some new semi-quantitative 
analysis. Our intent is to keep work analytically tractable at this stage. 



We set fcc = A = in ()2.12p and ()2.13|) . In the Einstein-deSitter model, the critical 
density can be written exphcitly as 

and the relation + = 1 holds. We also introduce quantities resembling the comoving 
Jeans wavenumbers for each component taken separately: 

= 11 = '^. (2.18) 

The difference with the true comoving Jeans wavenumber for a one-component fluid is the 
inclusion of the total mass density po, rather than just the mass density of the component 
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The effects of the various physical processes are now clearly evident. The expansion of the 
Universe produces a damping term 35-/(2a), causing the solutions to be in power law form 
rather than exponential. The relation of the mode wavenumber to the Jeans wavenumber is 
expressed as a ratio, transparently showing in which region of physical scales the mode lies. 
This ratio can be compared to the fractions and e^i to decide whether gravity or pressure 
dominates the dynamics. The true Jeans instability scale for a two-component medium is 
not given by either ks or ku, but by a combination of the two, as demonstrated in ^ [ij. 
This scale will be introduced in due course. 

We finally perform a couple more manipulations, to cast the equations in their simplest 
form. We define the dimensionless parameters 

Kb = 7^, Kd = (2.21) 

ks kf) 

Then Ki < 1 corresponds to the Jeans unstable region in the one-component analog of the 
equations, and i^'j > 1 to the acoustic region. We also make the variable transformation 
X = a^^/^. This gives the final form of the system of differential equations to be studied in 
the ensuing sections: 

S's + 6 (i^l - ^) 6b - ^Sn = 0, (2.22) 
6', + 6[Kl-'-^yn-'-^6s = 0. (2.23) 
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A prime now denotes differentiation with respect to x- These equations bear a strong resem- 
blance to the equations of an electron-proton cosmological plasma studied in [3] [equations 
(4.8) and (4.9) of that paper]. As is well known from the analogy between the simple Jeans 
instability and Langmuir modes, this resemblance is not surprising when the mathematical 
similarity between the electromagnetic and gravitational forces is considered. The tech- 
niques employed in 3] will be useful in our current analysis. In this paper we will employ 
the Frobenius method in obtaining long wavelength solutions. In a related paper some 
general WKB techniques are developed further than previously. Apart from facilitating 
some short wavelength solutions to the current problem, these techniques will also indicate 
further results possible in cosmological plasma physics. 

Before we proceed to a general analysis of (j2.22|) and (j2.23j) . we wish to digress to the 
simpler case of a one-component system. Surprisingly, we will derive some apparently new 
results, which provide a conceptually useful introduction to the ensuing analysis. 



III. THE ONE-COMPONENT EQUATION REVISITED 

The Einstein-deSitter one-component equation for a baryonic or dark matter fluid in 
the post-recombination era is a canonical example studied in all textbooks for linearized 
cosmological perturbation theory. It gives the familiar Jeans unstable power law solutions 
5 oc t^/'^, in the limit of large scales, and acoustic oscillations in the limit of small scales. 
Despite this, we have not found the full exact solutions completely displayed and analyzed 
in any textbooks or review articles in the current literature. Although a full analysis will not 
bring any startling new physical revelations, the mathematical techniques required are of 
some interest in their relation to the physics, and as an introduction to the more complicated 
analysis we will require later. This section may be seen as a useful orientation to the further 
work carried out in the bulk of this paper. 

We begin with the one-component version of ()2.19|) . i.e with = 1 and e/j = [or vice 
versa for ()2.20j) ]. Analogous to the definitions of Kb and J^d? we define the one-component 
comoving Jeans ratio for the fiuid, which has comoving Jeans wavenumber kj, as Kj = k/kj. 
The one-component density perturbation equation then becomes 
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The solution of this equation is a Bessel function of order 5/2. A Bessel function of half 
odd-integer order can be recast in terms of a spherical Bessel function. To begin with, we 
will choose the spherical Bessel functions of the first and second kind, jy and yy respectively. 
The solution may be rewritten as 



6{a) 



cia ^/^j2 




Kj 1 + C2a-'/^y2 




(3.2) 



with arbitrary constants of integration ci and C2- For this case the fortuitous circumstance 
arises that the solutions may be represented in terms of elementary trigonometric functions 



(see e.g. 2^). The solutions as shown are exact mathematical representations, containing 
all the information of the modes over all scales. As is usually the case with such solutions, a 
simple inspection does not reveal all the physical properties of the modes in an obvious man- 
ner. For example, it is a little difficult to interpret the time dependence of the modes through 
the argument of the Bessel functions \f^JaKj. We will require various approximations and 
numerical plotting to extract more physical meaning out of the solutions. 

To begin with, we seek to place the solutions into a canonical form, for easy comparison 
with other examples. The most useful such form comprises, to leading order, a product of 
a power law time factor and complex exponential factor. This approach was adopted in 



the studies of cosmological plasmas 
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18| for one- and two-component systems. As 



mentioned previously, due to the similarity between the gravitational and electromagnetic 
forces, the corresponding modes display many similarities. 

The most useful Bessel function solutions for our purposes are the Hankel functions, due 
to the fact that their leading order terms contain complex exponentials. We use the spherical 
Hankel functions /ig^"* and /i2^\ given by the expressions 



h 



h 



(2), 



1 

- exp 

z 

1 

- exp 

z 



37r\ 
_ 37r 



i-l-i 

z^ iz 



. 3 3 
1-- + - 

z^ iz 



(3.3) 

(3.4) 



We may make the comparison here to plasma results, where analogous series were obtained 
for large z. Contrary to here, where the series has a finite number of terms, the series for 
plasma modes were only asymptotic. 

We write down the explicit one-component solution via Hankel functions as 



6[a) 



1 + 



2K| 



+ 



1/2 



exp < ±i 



__ + arctan ^^j^ 



1/2 



1 



(3.5) 
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It must be stressed that unlike the plasma solutions, this is an exact result. The modulus of 
the solution grows with respect to time to leading order as 5 oc a if Kj <^ 1, or else if Kj ^ 1 
the modulus is approximately constant, with a first order time correction proportional to a. 
There is also a complex exponential portion to the solution, which usually gives a dispersion 
relation. The dispersion relation may be extracted by differentiating the phase with respect 
to t. This follows from the general fact that given an observed frequency u, a solution of 
the form 



S oc exp 



±i j uj{t)dt 



(3.6) 



is expected. This is assuming, of course, that the solution oscillates — if not, some other form 
of real valued solution must be available. Using the matter dominated time dependence of 
a, namely 



a = ( - ) , (3.7) 



where ti is an arbitrary constant, we find the frequency to be 



1 + \Kfa + \K^^a^ 



The result has been expanded for Kj > 1, as we may suspect that due to the Jeans instability, 
acoustic waves only exist in this region, and thus we can only attach physical meaning to lj 
for Kj > 1. This assertion will be derived rigorously in what ensues. 

The result of ()3.8|) may come as a surprise. How does it relate to the well-known Jeans 
dispersion relation derived for a static spacetime 

In a cosmological setting, we may expect the dispersion relation to follow a similar form, 
with appropriate time factors included. For plasma modes, it was demonstrated in that 
the dispersion relations could be written down to leading order in exactly the same form 
as their static spacetime counterparts in terms of physical (non-barred) variables, and then 
converted to comoving variables by inserting the correct time factors. Thus we may expect 
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at least in the form of a binomial expansion, namely 



UJ ~ 




2/n 



1 a 



1 




1 

16^ 



) 



(3.11) 



This form for uj may also be expected to contain some other time dependent terms, as 
was demonstrated for a number of plasma modes. Comparing the expansions in ()3.8|) and 
(13. mi . we see in fact that they only agree to first order. This still indicates some form of 
Jeans instability, but the dispersion relations are quite different. This difference in behavior 
between the linearized gravitational modes and plasma modes may be attributed to the 
special role the density plays in the gravitational perturbation equations. Eq. p.ip contains 
only one free parameter, the Jeans ratio Kj^ whereas the plasma equations contain both 
the sound velocity and plasma frequency, which cannot be reduced to one parameter. This 
implies that the relation between the gravitational source and the Friedmann equation, 
which fixes the background spacetime, means that the same form for the dispersion relation 
as found in static spacetime need not necessarily be expected in the expanding Einstein- 
deSitter model. 

We have a general solution in terms of a modulus and complex exponential, which is 
exact and thus contains all the information of the problem. How do we infer the usual Jeans 
instability behavior from this? Let us examine plots of the solutions to gain a pictorial 
idea of what is happening. In Figs. Q and El we see the transition from acoustic oscillations 
to growing and decaying modes as Kj is decreased — we are examining ever larger scales, 
and passing through the instability. To all fit on the same set of axes, the plots have 
been approximately normalized. The dependent variable a is rather arbitrary [as can be 
deduced from ()3.7j) ]. An appropriate starting time ti may be chosen to normalize a to 1 at 
the beginning of the chosen epoch of evolution, and the solutions may then be propagated 
forward in time. The tendency for the period of the acoustic oscillations to grow longer in 
time is evident from the plots, and the approximate constancy of the amplitude predicted 
previously from ()3.5|) is evident. In the extreme case, the oscillation period becomes so long 
that a perturbation cannot complete one full oscillation, and then the instability arises. This 
behavior can clearly be seen in Figs. Hand|2lfor the Kj = 2.0, 8.0 plots. 

We now perform some approximations to make contact with some better known results of 
the one-component problem. Let us begin with a small Kj/a expansion. For the spherical 
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Hankel solutions ()3.5p . we find 
6 



2Kj 



exp 



al/2 



5 a2 ^ ' ' 



(3.12) 



The above expansion explains what happens to acoustic oscillations when Kj < 1. In 
this region, the leading order factor Kja~^^'^ in the exponential must always lie between 1 
and numerically, and decreases with increasing time. This is because a > 1 and increases 
monotonically for all time. Thus the solution lies within one period of oscillation for all time, 
and only the growing or decaying modes may be observed. When Kj becomes larger than 1, 
more than one period of oscillation may be spanned by the Kja~^^'^ factor, and the solution 
will begin to develop acoustic waves. The expansion ()3.12|) shows 5 oc a, which only gives 
the familiar growing mode, discussed in all texts. The decaying mode has not been found 
in the current analysis, because spherical Hankel functions have been chosen to represent 
the Bessel function solutions. The spherical Hankel functions are linear combinations of 
the original j2 and 1/2 solutions, which contain both modes. The decaying mode has been 
"asymptotically swamped" by the growing mode in this linear combination. The decaying 
mode may be liberated by a direct small variable expansion of fj3.2j) for each spherical Bessel 
function. We find the j2 component gives the decaying mode 



5 ~ 



7 a V 



(3.13) 



and the j/2 component gives the same growing mode as ()3.12j) in the slightly different form 



S 



1 + ^ + { ^ 
a \ 



(3.14) 



without the exponential phase factor. These solutions clearly correspond to the usual text- 
book modes found when pressure is ignored. They have included the pressure corrections, 
given as a series in the Jeans ratio, with matching time factors to all orders in the expansion. 

We now turn to the large parameter expansion, from which we expect to liberate the 
acoustic oscillations. Eq. ()3.5|) is in fact already in the form of a large parameter expansion, 
only that it has a finite number of terms, and is consequently exact. The phase of the 
exponential does not seem to show the structure of the familiar Jeans dispersion relation. 
This was discussed above, where it was indicated that this is not necessarily to be expected. 
It is possible to see why this is so in a lucid fashion by applying the WKB method to the 
original equation ()3.1|) . Through this method we will derive a dispersion relation displaying 
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very similar characteristics to the famihar textbook one, that shows the presence of the Jeans 
instabihty. A more comphcated WKB approximation scheme was developed in to deal 
with plasma modes of a more intricate form, involving larger numbers of coupled equations. 
This method is used in [19] to handle the two-component cosmological density perturbation 
equations in the short wavelength approximation. For the present simple second order 
equation, the standard textbook approach suffices (for a good explanation of WKB methods, 
see I29I). 

To see the physics most clearly, we transform ()3.1|) to depend on t. Using the variable 
transformation given by (|3.7|1 . the equation 

results. The usual Jeans dispersion relation can be directly seen in this equation. Consider 
the factor 2/(3t^), arising from the gravitational source term. This term may be transformed 
to explicitly see the source parameters emerge. We consider the physical (time dependent) 
value of the total energy density: 

4.Gp„ = 4.Gp.= ^ = |,. (3.16) 



and the relation 



Kj = lvrkHl (3.17) 
Then the Jeans dispersion relation can be directly seen in ()3.15|) : 

vrk'(jy' -^ = v!k^-A7rGp,. (3.18) 

In the static spacetime case the physical variables would of course not depend on time, and 
the first derivative term 45/ (3t) in ()3.15|) would not exist. This leads to the exact exponential 
solutions and the familiar dispersion relation given by ()3.18p . as originally found by Jeans. 

The most straightforward way to effect a WKB approximation in the present situation 
is to remove the first derivative from the equation. One way to do this is by the variable 
change x = a^^^'^- We then find the equation 

'l^^^ + (6Kj-^)s = 0. (3.19) 

Now we suggestively define 



Q X 1/2 

u{x) = { QKj - . (3.20) 
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Applying the WKB approximation to ()3.19|) gives the leading order solution 



6{x) ~ ^ 



~-l/2 



exp 



(3.21) 



A dispersion relation has indeed been derived, and is given by oj as defined in (j3.2Up . If we 
consider 



1/2 

- ^\ dt 



8/3 3^2 



-4/3 



3tf \t 



2/3 



1 1/2 



dt. 



(3.22) 



then a physical frequency uj(t) can be identified by use of ()3.18|1 and the relation tf 
l/(67rG'po). Thus 



'^{x)dx = 



i{t)dt = j {vlk'^ - AnGpoY^^dt, 



(3.23) 



and the physical connection has been made. Another point to note is that the amplitude 
a}(x)^^/^ is time independent to leading order (because Kj ^ 1). Thus the amplitude is 
approximately constant, as noted previously. 

To make a direct comparison of ()3.2H) and ()3.5|) . we now evaluate the integral in the 
phase of the WKB solutions. A change of integration variable from x back to a puts the 



integral into a form which has been tabulated 

1/2 



and we find 



i^{x)dx 



QKj 



6 



arcsin ( 1 



2a 



(3.24) 



The solution generated by the WKB method may be compared to the exact one given by 
()3.5p . The amplitudes and phases need to be expanded for large Kf. 



WKB phase: 



-6 



1/2 



arcsin ( 1 — 



2a 



Kj 



f.V2 



TT a 



1/2 



exact phase: 



+ arctan 



4 Kj 

ja 



(3.25) 



\2K\ 



1 



v/6 



1 + 



a 



2Kj 



+ 4, 



WKB amplitude: {<oK] - 6a)-^/^ 



gl/4;^l/2 



1 + 



exact amplitude: 



1 + 



a a 
2Kj^4Kj 



1/2 r 



1 + 



4Kj 



(3.26) 
(3.27) 

(3.28) 
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It can be seen that the two solutions only agree to leading order (modulo time indepen- 
dent constant factors). Given that the WKB method only gives a leading order solution 
to the problem, no more can be expected. This discussion has highlighted the difference 
between the expected dispersion relations of static spacetime to those derived in an expand- 
ing universe context. WKB can reproduce the same form as the static spacetime dispersion 
relations, but this may only agree to leading order to the true dispersion relation found in 
an expanding universe scenario. 



IV. IMPROVEMENTS ON PREVIOUS CDM PERTURBATION RESULTS 

We now return to the two-component equations and consider the case of CDM pertur- 
bations characterized by strictly zero temperature. If the velocity dispersion is considered 
to be an adiabatic sound velocity as given by (j2.8|) . then Tj = corresponds to Kj:, = in 
()2.2H|) . Such an approximation facilitates an exact analytic solution to the problem, which 
is otherwise impossible. The resulting system of equations was one of the main cases in- 



vestigated in the general analysis of In this section we will point out what appears to 
be an error in the analysis of that paper, which leads to some markedly different solutions, 
derived in what ensues. 

Before we proceed with this, it is pertinent to point out a general problem with taking 
Ko = in Eqs. ()2.22j) . ()2.2H|1 . For simplicity, it is possible to neglect spacetime expansion, as 
the qualitative behavior will be the same. Thus the static spacetime cosmological equations 
studied in ~\ are sufficient for this discussion. These equations are also examined in jioj ]. 
and using the notation employed there we have 

5d + {vie - Wd)5d - Wb5b = 0, (4.1) 
'^B + {vie - Wb)6b - Wd6d = 0, (4.2) 

with Wi = inGpi. This system of equations can be reduced to a first order linear autonomous 
dynamical system describing a state vector 

X = {xi, X2, X3, x^)^ = {6d, Sd, Sb, Sb)'^, (4.3) 

where T denotes the transpose of a vector. The dynamical system has a /c-dependent critical 
point found by solving the equation x = 0. This critical point happens to give the Jeans 
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wavenumber of the mixture, and is given by 



e = ki, ^ki + ki = '^ + '^. (4.4) 

Note the shght difference in this definition of ki in comparison to ki defined in fl2.18|) . It is 
aheady evident that a problem arises if we take f^) — > in ()4.4|) . and this will be physically 
elucidated by studying the four independent modes of the system, given by the eigenvalues 
of the dynamical system. 

The eigenvalues give the structure of the modes (see [3] or [3] for the details). They are 
found to be of the general form 



For the case of CDM currently under consideration, where vd = 0, the fc-dependent functions 
/ and g are given by 

f{k) = Wb + Wd- k\l, (4.6) 
9{k) = ewlvl. (4.7) 

It was previously found that the eigenvalues Ai and A2 described the Jeans unstable modes, 
whereas A3 and A4 described acoustic oscillations at all wavenumbers. An examination of Ai 
in the current context will show this not to be the case for CDM. If we define (analogously 
to the one-component scenario) 

Wb + Wd ^ ' 



Ai may be written as follows: 

1/2 

V2 



Ai = ^{Wb + WoY" {l-K]+[{l- K-f + Ae^K^j] . (4.9) 



This does not equal zero for Kj = 1, and it is straightforward to show that it has no zeros 
for all 7^ 0. Thus a CDM perturbation would collapse for all scales, clearly a physical 
impossibility. 

If we examine Eqs. (j4.1|) . (j4.2p . it can be seen that with the removal of the pressure term 
V£,k'^, there is no mechanism to counter the remaining gravitational source terms —Wo^d 
and —WbSb, whose sign indicate an attractive forcing, initiating gravitational collapse. This 
is the case no matter how small the fraction of dark matter compared to baryons, thus no 
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amount of baryonic pressure support can prevent a collapse at any scale. This physically 
absurd situation is the root of the problem of taking vd —>■ in the fluid models. A physically 
correct equation must include some sort of velocity dispersion term, even if the matter is 
totally collisionless and an adiabatic speed of sound cannot be defined. 

With these thoughts in mind, we must view the current section as more of a mathematical 
digression, than a physically realistic model. It nevertheless serves a purpose. We can 



make contact with the work of and uncover some interesting mathematical properties 
associated with these cosmological perturbation equations in general. We find mathematical 
subtleties overlooked in 0| , which also indicate the nature of the general solutions to follow 
in the next section. They allow an interesting comparison with cosmological plasma modes 
discussed in [18], where the limit — > is valid. The special nature of gravity, and the 
extra complications it entails are revealed by this comparison. 

We now proceed to obtain a solution of Eqs. ()2.22p . ()2.23p (with Kd = 0) in the long 
wavelength (small k) limit. A short wavelen gth solution is trivially obtained by setting 
Kd = everywhere in the results presented in . To align ourselves with earlier notation, 
and to stress the fact that there is only one Jeans related scale now occurring, we will 
rename Kb to Kj. Rather than directly reducing Eqs. ()2.22|) . (12.231) into a single equation, 
we attempt a solution by the Frobenius method. This is useful as a precursor to the general 
solution derived in the following section in a similar manner. 

To begin with, we assume an arbitrary expansion of the solutions in the form 

oo 

Sb{p,x) = x''5Za„x", (4.10) 

n=0 

oo 

n=0 

Here p is an arbitrary exponent to be determined, and a„ and 6„ are series coefficients 
also to be determined by the Frobenius method. We substitute these series into fl2.22|) and 
fl2.23|) to obtain a set of algebraic relations between the undetermined coefficients. Then all 
coefficients of like power of x collected and equated to zero. 

Arising out of this procedure are a pair of indicial equations for p, with one of qq or bo 
remaining an arbitrary constant: 

p{p - l)ao - Gesao - Geobo = 0, (4.12) 
p{p - l)bo - Geobo - 6eBao = 0. (4-13) 
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A set of recursion relations also arise for all higher order coefficients: 

6Kjan +[{p + n + 2){p + n+l) - 6eB]a„+2 - QeDbn+2 = 0, (4.14) 

{p + n + 2){p + n + l) - 6eD]6„+2 - 6eBa„+2 = 0. (4.15) 

In retaining the arbitrary constant Qq Bq, it is assumed that all odd indexed terms vanish 
from propagation of the initial values ai = bi = through the recursion relations. 
If we solve (|4.12|) and (j4.13|) for p, we find four possible values: 

p = 0, 1, 3, -2. (4.16) 

A comparison with the analogous case for cosmological plasma modes immediately 
shows a difference in the nature of the exponents. In the present case the exponents are 
exactly determined integers, whereas for plasmas the exponents depended on the plasma 
frequency. When the one-component solutions were discussed in the previous section, an 
analogous difference was observed between the spherical Bessel function solutions of the 
gravitational perturbation modes, and the general order Bessel function solutions of the 
plasma modes. Whereas in the one-component study the solutions were simplified by this 
property of gravity, in the present case they are in fact complicated. The plasma solutions 
were representable in terms of 2-F3 generalized hypergeometric functions, but a similar rep- 
resentation is not well-defined here. This is because the exponents differ by integers — a fact 
which necessitates a modification of the basic Frobenius method. This modification is borne 
out in the solutions by the fact that parameters appearing in the denominator of general- 
ized hypergeometric function expansions cannot differ by integers. In such situations the 
generalized hypergeometric functions are not definable, and one must resort to classifying 
solutions of the equation by Meijer G-functions. 

To apply the Frobenius method to indices differing by integers, the recursion relations 
must first be solved for general p. This is achieved by employing ()4.12|) . ()4.13|) . and writing 
the system of differential equations as 

p{p - l)ao - QeBtto - Qeobo 
p{p - 1)60 - 6ez)6o - 6eBao 



L 


^b{p,x) 






Sd{p,x) 




where the operator L 


is defined 


by 



X 



p-2 



(4.17) 



6e 



D 



6eB 
' X' 



X 

92 6e 



D 



X' 



(4.18) 



21 



At present there is no relation between oq and bo, but the recursion relations will provide 
one. A lengthy algebraic exercise is needed to solve the recursion relations. The result is 

~ ^'^D + |)n(|p + ^f^D + |)n 



J2n 



(|p - \)n{\p + \)n{\p + l)n(iP + 2)„ 



(-1^3)' 



p{p - 1) - 6eD C2P - + + IUIp + 2) 

In the above, the parameter has been introduced as a shorthand: 



(4.19) 



(-1^5)". (4.20) 



(4.21) 



and the notation ()„ is the Pochhammer symbol. It is clear that the coefficients as derived 
will not exist for p = 1, —2. This is the basis for requiring a modification to the straightfor- 
ward method of substituting in the four calculated values of p ()4.16|) into ()4.19|) and ()4.2U|) 
to generate four independent solutions. 

To proceed, we extract from (j4.19|) and (j4.2m) the relation 



which allows (14.1711 to be rewritten as 



Snip, X) 
Sd{p,x) 



ao 



p{p- 1) -Q(^d'' 



1 p(p-l)(p + 2)(p-3) 
es p{p - 1) - 6eD 




(4.22) 



X 



p-2 



(4.23) 



For p = 0, 3, direct substitution is permissible to obtain two independent solutions, which 
turn out to be generalized hypergeometric functions, analogously to the plasma results. 
For p = 1, —2, we take advantage of the fact that ao is arbitrary, and set it respectively equal 
to p — 1 and p + 2. After differentiation with respect to p and evaluation at the respective 
points p = 1 and p = —2, we obtain the result 

d 



-^Bip^X) 



a 

L^^..(P,X) 



0. 



(4.24) 



p=l,-2 



This gives two more solutions. 

The above discussion has defined an algorithm for finding small k expansions of the 
solutions of the cosmological density perturbation equations. This algorithm is readily im- 
plemented into a symbolic manipulation computer code. We will write such a code to find 
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the solutions to the general equations in the next section. The solutions for the CDM case 
under consideration may now be written down. The need to differentiate with respect to p 
for the p = 1, —2 cases result in solutions involving digamma functions ip and logarithmic 
;erms — quite a complication to the 2-^3 functions derived for the plasma perturbations in 
isl ] . To display the modes in their simplest form, linear combinations of the modes derived 
directly by the above algorithm need to be taken, and the use of various mathematical iden- 
tities involving digamma functions employed. The final set of four CDM modes are given as 
follows: 

" (4).4)»(2)»"! *" ' 

= - 2-Fs (f - \vD, I + \vd; -i i 2; -|A'3x') (4.25) 



^ T? f 1 1,, 111,,. 11 

en 



2F3 (-i - '^un, -i + lun; -i, i, 2; -iKjx') (4.26) 



^ IsJnU''"'^"' ~ -'-''•'^• 
X [V' (f - I'^D + n) + (I + |z/B + n) - - ^(n + 1) 

-ij{n + I) -ij{n + I) +\ogx'] (4.27) 

X (i - I'^i^ + n) + (i + |z/B + n) - V'(n) - ^(n + 1) 

- V^(n + I) - ^(n + I) + log x'] (4.28) 

'5b3(x) = Csx\F, (f - I + I^^d; I, i, 2; -fi^5x') (4-29) 

5d3(x) = cax' 2^3 (f - I^^D, f + I^^d; |, i, 2; -f i^^x^) (4.30) 

X [V- (I - I'^D + r^) + ^ (I + I^^D + n) - V^(n + 2) - V^(n + 1) 

-ij{n - I) -^{n + I) +\ogx'] (4.31) 

5.4(x) = C4 iiKjxr' - ^c.es ± (-1^3^^)" 

X [V' (-i -^UD + n) +ij (-i + iz/i? + n) - ilj{n + 2) - ^(ri + 1) 
-V(n-i)-V(n + i)+logx'] . (4.32) 
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The same solutions may be obtained by considering Meijer G-function solutions to the 
original differential equations. The procedure for determining solutions to generalized 
hypergeometric-like equations which contain parameters differing by integers is discussed 
in detail by Luke l^l, pp. 138-143. The study outlined a method for developing Meijer 
G-function solutions from the equations, which involved the differentiation of generalized 
hypergeometric functions with respect to their parameters — a procedure analogous to the 
differentiation of p indices in the above. The analysis involved is lengthy and tedious, but 
leads to the solutions obtained above. We refrain from a physical interpretation of the above 
gravitational modes for now, and take that up in the next section when we discuss the more 
general solutions. The one obvious difference will be the fact that no Jeans instability is 
apparent for any of the above modes, with two modes always collapsing and two modes 
always acoustic, whereas for the general modes a Jeans instability will be apparent. 

We now compare the solutions obtained to those of In these studies general 

solutions were written as 

01 = = ciGi + C2G2 + C3G3 + C4G4. (4.33) 

Here q are constants, and the functions Gi denote Meijer G-functions 



Gh — G2 4 I X 



1, 2, 3, 4. (4.34) 



The notation # signifies that is to be omitted in its usual place. The G-function is defined 
so that < m < 4 and < n < 2. The parameters a*,h* depend on e^, ^d, the adiabatic 
index 7^, and the exponent of t in the Hubble expansion parameter (which we pointed 
out earlier must be equal to | for the equations as formulated to be physically correct — even 
though 0, 11 1 used a greater range of values). The time parameter x = ^KjX^ in our 



notation. A particular set of solutions is given by m and n being given specific values. In 
general, m = 1, n = 2 will give such a set of solutions for small x in the above example. 
The G-functions can then normally be expressed in terms of 2-^3 functions (for example the 
plasma solutions), but in the particular case under discussion, since some of the b* differ by 
integers, this is not possible. It is this point that Haubold and Mathai missed in [ll|. 

Let us get down to specifics to illustrate the point. Under the general classification 
scheme, the CDM case under consideration corresponds to 77 = |, 7^ = |, in the notation of 
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IjJ. This implies that the parameters take the following values: 

al, a*2 = -1± iz/D, 

To evaluate the G-functions in this special case, where some parameters differ by integers, 



the integral representation of the G-functions was considered in 



1 f T(^ + s)T(-a; - s)T(-a%- s) 



Using the fact that 



it was found that 



n-\-s) n{-s\ 



(4.36) 



c - ^5/4 r(-at + !)r(-4 + |) p ( . , 5 5. 5 7. \ 

Gi - -X — r(|)r(2)r(i) + i' + 4, 2' 2' 2' "V ' ^^"^^^ 

which agrees with in ()4.29|) . Performing the same analysis as above for however, it 
is found that Gi and G3 are exactly the same function. It appears as though this was never 
checked in The same applies to G2 and G4: 



r - r - .-i/4 rMi-i)r(-a2-i) 
''^-''^-'^ r(i)r(2)r(-i) 

X 2F3 [-al - ^, -a; - ^; ^, 2, -x ) , (4.38) 



which agrees with in ()4.25|) . Thus the solutions degenerate by employing this method, 
and two of the solutions, namely 6b2 and 5^4 which contain logarithmic terms, are totally 
missed. The correct way to evaluate the G-functions when some parameters differ by integers 
i. .iven ,„ Lu.e H PP-MS-MT. A ca.efu, .M. of .M. .... co.p„ca.ed p.o.d.e w,U 
show that the logarithmic solutions 6b2 and 6b4 are found in this way. 

It is interesting to note that a linear combination of the solutions 632 and 6ba are in 



fact found in ll| in the large x limit. In this case the function 6*24 ^as evaluated from 
the contour integral representation to achieve some analogous series to ()4.27|) and ()4.31|] . 
Although the solutions ()4.25|1 - ()4.32|1 are valid for all x, such a representation does not seem 
very useful for x large, as the solutions are comprised of an infinite ascending series in x. 
Thus very many terms would be required to represent the solutions accurately in this limit 
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through 5b2 and Sba- In (l9| we develop a much better method for evaluating the solutions 
for large x using a WKB approximation scheme. 

We have now thoroughly investigated the CDM two-component model using the K^j = 
limit, and tidied up previous work in this area. As discussed earlier, although this work 
is of dubious physical relevance, it has been an interesting mathematical investigation, and 
has allowed comparisons with the previous work in cosmological perturbation theory and 
cosmological plasma physics. We now turn to the more general two-component model, which 
has a firmer physical basis. 



V. THE GENERAL TWO-COMPONENT SOLUTIONS 

We finally investigate the most general set of equations, those with both Kb and Kd 
nonzero. The equations thus posed can model both CDM and HDM, though the discussion 
initiated earlier about the lack of a free streaming damping term for neutrino-like HDM 
should be heeded. 

The algorithm needed to solve the equations was developed in the previous section for 
Kd = 0, and can be used here. Thus if ()4.10|) and ()4.11|) are substituted into ()2.22|) 
and ()2.23|1 . we obtain the same indicial equations ()4.12|) and ()4.13|1 as previously, and the 
following set of coupled recursion relations: 

{p + n + 2){p + n + l){p + n + 4:){p + n - l)a„+2 
+ GKKp + n + |-z/z?)(p + n + | + UD)an + SGeoKlbn = 0, (5.1) 

{p + n + 2){p + n + l){p + n + 4:){p + n - l)6„+2 
+ 6Kl{p + n + |-z/B)(p + n + | + i/b)^™ + 36esi^'|a„ = 0. (5.2) 

Here z/^ = | -|- Ge^ in analogy with the previous definition of I'd- A closed solution for this 
set of recursion relations cannot be obtained for all n, so that solutions to ()2.22|) and ()2.23|) 
can only be generated iteratively, to whatever order desired. Lacking the ability to generate 
an infinite series representation for the solutions means that they cannot be classified by 
known analytic functions. To handle the complicated algebra involved in finding successive 
terms iteratively, we have developed a symbolic computation code using the functional 
programming language Mathematica. The algorithm described in the previous section can 
be used to generate a solution up to a certain power in x- 
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The coefficients increase in complexity very quickly for increasing n. Although the code 
can generate solutions up to arbitrary order, we find it sufficient to present only the first 
two orders for each solution here. We express the results in terms of the original parameters 
Kb and i^D, rather than in terms of Vb and as no simplification is gained in using the 
latter. The solutions, corresponding to p = 0, 1, 3, —2 respectively, are: 



^Di(x) 



1 + -{Kl - SenKl - 3eBKl)x' + O(x') 



CD 



1 + -{Kl - SeoKl - 3eBKl)x' + 0{x' 



(5.3) 
(5.4) 



1 
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X+-{QKl- SleoKl - ^UbKI + h'-^Kl ) 



CD 



+ ^-eB{Kl- Kl)x'\ogx + 0{x') 



8d2{x) = 



CD 

6 



X - - {SlenKl + 3UbKI - Kl) 



-en{Kl-Kl)x'logx + 0{x') 



(5.5) 



(5.6) 



^Bsix) 


= X' + ^{^^dK'b - SeoKl - lOKl 


)x' + o{x'), 


(5.7) 


Sosix) 


= x' + ^{S^bKI - 3eBKl - lOKl 


)x' + o{x'), 


(5.8) 


^Biix) 


= X-' + (kI - 5'-^Kl + 5enKl - 

V CD 

+ QeDiKl-Kl)\ogx + Oix'), 




(5.9) 


SoAix) 


= x'^ + {QKl - 5eBKl + beBKl) 
+ 6eD{Kl-Kl)logx + Oix')- 




(5.10) 



It remains now to make a physical interpretation of these solutions. This is most usefully 
achieved by constructing a table comparing different properties of the modes, and comparing 
the solutions to the corresponding cosmological plasma modes. The results are summarized 
in Table lU which uses notation defined in jl^ . Both the Kb = modes of the previous 
section and the modes of this current section are included in each category under the table. 
The plasma modes yi . . .y^ of [l^ Eq.(4.16) correspond to the Kb = modes ()4.25|) - ()4.32|) . 
The current gravitational modes ()5.3|1 - ()5.1()|1 (that is (5i . . . ^4) correspond to the more general 
expansions Eqs.(4.20), (4.21), (4.29)-(4.32) of yj]. 

In Table E the power of X gives the corresponding exponent p. The parameter u = 
\ — Pf — Pe found in the plasma modes depends on the plasma frequencies of the electron 
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TABLE L A comparison of gravitational and plasma linear perturbation modes 

and ion components. For the gravitational modes, the signs of Pi and Pe must be reversed 
since (as opposed to the electromagnetic force) gravity is always attractive. In addition 
P^ + P^ corresponds to €b + €d — i (once again the special nature of gravity in cosmology is 
apparent). This implies that the general parameter u in the plasma modes should be replaced 
with I for the gravitational modes — the reason why many of the parameters in the 2-^3 and 
Meijer G-functions were pure rational numbers not depending on physical constants. Notice 
that the ratio of the amplitudes of baryonic/dark matter modes and electron/ion modes 
differ. This is because the couplings in the differential equations are different. For the 
gravitational modes the couplings involve terms such as sdSd and cbSb, whereas for the 
plasma modes the couplings involve terms such as Pfuei and P^un- 

We have indicated the corresponding collapsing and acoustic modes for the gravitational 
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density perturbations. It is difficult to show this rigorously for the series solutions as pre- 
sented. We can make comparisons to the one-component results, and identify the leading 
order powers of the expansion parameter a. This yields the classification as stated. We can 
also make an analogy to the ion-sound modes of plasma physics, which are of a similar na- 
ture to acoustic oscillations. They show a collective behavior of both components oscillating 
approximately in phase. 

We are finally left with the question of how the Jeans scale enters into the solutions. In 
Section 4 the mixture wavenumber fc^, Eq. ()4.4j) was briefly introduced as being the only 
physically meaningful scale for instabilities in a two-component fluid. To make this quantity 
dimensionless, it would be appropriate to make the definition 

This quantity is only of relevance to a static spacetime scenario. To place it in the context 
of the expanding Universe, the substitutions 



Wb- 


6eB 




Wd- 




vie QKl 




X 



are required. Then Km takes on the revised definition 

,2 ^,2 



~ I TZI I I T^l ~ 9" (5.12) 

We have introduced the quantity Xc{k)-, which can be thought of as a critical time. For 
X > Xc, Km > 1 and acoustic oscillations would only be expected to exist for all modes. For 
X < Xc, Km < 1 and two of the modes become unstable and undergo gravitational collapse. 
Since the precise magnitude of the scale factor a is not determined by cosmology [see Eq. 
fl3.7|) . which contains an arbitrary initial time tj], we may arbitrarily assign an initial time 
= 1, so that Xo = 1 and decreases with increasing time. Then we may interpret the 
Jeans instability in two ways by considering the critical time Xc- Initially we may study 
all fc-dependent modes at a particular time x, where a subset will be unstable for values of 
k for which Xc{k) > X (we stress that Xc is a function of k). We may then consider what 
occurs as the modes evolve through time from this particular instant. The critical time Xc 
is fixed for any particular mode, so that a subset of modes that were originally acoustic 
will become unstable &s x ^ xt (those modes corresponding to the solutions ^3 and 64 in 
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Table I). Consequently more and more modes pass through the instability as the Universe 
evolves. The physical wavenumber k is of course dependent on time, thus the dependence of 
the instability on a time Xc shows the inextricable link between the wavenumber and time. 

It is illuminating at this stage to refer back to the one-component modes discussed in 
Section 3. For the one-component case, solutions were found in terms of the combination 
Kja'^^"^. With the identification x = ar^l"^ ^ the quantity Xc is seen to be the two-component 
analog of Kj. 

It would be useful to convert the expansions ()5.3|) - ()5.10|) to depend on Km or Xc, to see 
how the Jeans scale enters. This is achieved by the following relations: 



1 



Xc 



(5.13) 



Here V = vb/vd is the ratio of sound velocities. We then find that the expansions are all 
in terms of increasing powers of x/Xc, with coefficients in terms of e^, cd and V: 
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X"' {1 + [-4eB + eoV^ + heseo + ^eliV"" - 1) 




+ Q^D (es - ££) + Ci)!/^ - -^j — \ogx + I — 



xi \xt 



) 



(5.20) 



^D4(x) 




(5.21) 



This is a convenient parameterization of the solutions. The scale of the modes are chosen by 
Xc, the nature of the matter involved is determined by V, and the proportions are determined 
by and eo- A complete solution to the problem has thus been achieved up to whatever 
order desired. 

VI. CONCLUSIONS AND FURTHER WORK 

A method for determining the small k solutions of a general two-component cosmological 
density perturbation model has been expounded in this paper. We have only displayed the 
solutions to first order, but it is possible derive them up to any order by the method in 
principle. We have explored the mathematical properties and peculiarities of density per- 
turbations influenced by gravitational interaction, particularly contrasting them to plasma 
modes, and correcting a number of previous misconceptions in the literature. The expanding 
Universe introduces new features not predictable from simple static spacetime considera- 
tions. In particular, totally new structures to the dispersion relations are found, even in the 
one-component example. We have shown how the mixture Jeans wavenumber enters the 
solutions, and clarified its role in an expanding universe context. 

More work is required to investigate the solutions around the critical scale defined by kM- 
Although the expansions as derived in this paper are applicable to this region, they are not 
particularly useful, as many terms in the equations need to be retained when the expansion 
parameter x/Xc is of 0(1). It is unclear how an analytical investigation of this region could 
proceed at present. We have performed some preliminary studies which involved producing 
a large number of terms in the expansions ()5.H|1 - ()5.1()|1 using the Mathematica program 
described, and then substituting in numerical values for the various physical parameters to 
obtain numerical coefficients with an ascending series in x- At present the plots of these 
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expansions over a range of values of x do not yield reliable results — it is possible that 
many more terms than are practically calculable will be required, and a very high order of 
numerical precision will have to be maintained. Other methods of analyzing the modes in 
this interesting region probably need to be investigated. 

Of ultimate interest is exploring how these type of modes contribute to the power spec- 
trum. More physical effects may need to be introduced, such as a cosmological constant, 
or the addition of more matter components. To determine the actual density contrast at a 
given scale the Fourier modes of the density contrast as derived in this paper would 
also need to be integrated over the whole range < A; < l/fc. It would be of considerable 
interest to compare the power spectra calculated by such a method with the well-known 
power spectra of the various cosmological models in existence today. 

In concluding, we remark that a similar analysis could be carried out in the post- 
recombination region 140 < z < 1150, where now the baryons follow the T Ija 
lationship. The differential equations in Section II will now be different, as will be their 
solutions; but we expect that the ensuing analysis would yield qualitatively similar results 
but quantitatively different scaling. This would be a useful future study. 
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FIG. 1: The transition of the decaying one-component modes as Kj varies through the Jeans 
instabihty 
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FIG. 2: The transition of the growing one-component modes as Kj varies through the Jeans 
instabihty 
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